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Abstract 

Exact solutions for nonlinear Arrhenius reaction-diffusion are con¬ 
structed in n dimensions. A single relationship between nonlinear dif- 
fusivity and the nonlinear reaction term leads to a nonclassical Lie sym¬ 
metry whose invariant solutions have a heat flux that is exponential in 
time (either growth or decay), and satisfying a linear Helmholtz equa¬ 
tion in space. This construction extends also to heterogeneous diffusion 
wherein the nonlinear diffusivity factorises to the product of a function of 
temperature and a function of position. Example solutions are given with 
applications to heat conduction in conjunction with either exothermic or 
endothermic reactions, and to soil-water flow in conjunction with water 
extraction by a web of plant roots. 
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1 Introduction 


Assuming that volumetric heat capacity is constant and neglecting reagent con¬ 
sumption, the temperature of a reactive mixture satisfies a nonlinear reaction- 
diffusion equation 


= V. [D(0)V0] + R{0), r e U C t e [0, ta), 


( 1 ) 


where U is compact and connected, V is the usual gradient operator, and t 2 > 0. 
R{9) is real-valued at any value of absolute temperature 9 G [0, oo). In order 
that heat conduction contributes to increasing entropy, D must be positive (e.g. 
m)- One of the most important forms for the reaction function is the Arrhenius 
reaction term 


r 9>0, 

\ 0 9 = 0 


with B a positive constant. This follows from the Boltzmann-Gibbs equilibrium 
canonical distribution governing the probability of a particle overcoming the ac¬ 
tivation energy barrier A of a reaction. The parameter B may then be identified 
with E/ks where ks is Boltzmann’s constant. The rate factor pCRo, p being 
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density and C being specific heat at constant volume, is the amount of heat 
energy released per unit volume per unit time at very high temperature. This is 
usually described by a weak power-law dependence on temperature, i?o = SqO^. 
For example, kinetic theory of a hard-sphere gas predicts m = 0.5 [5]. Since 
the temperature-dependence is dominated by the Arrhenius factor, most ex¬ 
tant models assume m = 0. This will be assumed here, but a straightforward 
generalisation of the following analysis can cover the case m 7 ^ 0 with pC also 
depending on 6. 

For first-order reactions, E is the energy of bond dissociation into reactive 
molecular components such as free radicals. For example, for the exothermic 
decomposition of diethyl peroxide, measured reaction rates agree well with the 
Arrhenius law over a broad range of temperatures [5]. In particular, as a com¬ 
bustible mixture is controlled by extracting heat through the boundaries of the 
container, the Arrhenius reaction form is expected to remain appropriate. 

A full Lie point symmetry classification of Q was made by Dorodnitsyn et al. 
[1]. Various combinations of power-laws and exponential functions for D{9) and 
R{d) lead to an expansion of the Lie group of point symmetry transformations of 
Q beyond the common Euclidean isometries in space, and translations in time. 
These in turn open possibilities of invariant solutions that may be obtained by 
reduction of variables. The invariant solutions indicate a wide range of possible 
dynamical behaviours, including stable similarity forms for the temperature 
with multi-peaks, extinguishing, and blow-up in finite or infinite time that can 
occur even in solutions with compact support [i]. 

The non-analytic expression exp(—B/0) is usually approximated by an ex¬ 
ponential function, in the Frank-Kamenetskii approximation [5], before it can 
lead to useful exact solutions of 0. Although unbounded reaction terms lead to 
insight on critical parameters for ignition, there is a significant difference in be¬ 
haviour when the reaction term is bounded, as in the standard Arrhenius model. 
The relationship between the full Arrhenius model and the Frank-Kamenetskii 
model can be more conveniently analysed by using the identity 


^-E/kB0 ^ ^-EtkBOo 


0 


1-f £0 


e 


E 


where 


0 = (0 - Oo) 


E 

ksdn 


is a rescaled temperature with the zero point shifted to some local value Oq (e.g. 
i)- Wake and Bazley [7] extended the lower bound 


^-E/kBS > ^-E/kB0o exp ( -^^ 

\ 1 -|- eOmax J 


beyond the maximum temperature to obtain close approximate critical values 
of parameters at ignition. Gustafson [5] constructed steady linear and radial 
solutions of the full nonlinear equation by “a shooting method combined with 
a Newtonian-Raphson technique and certain boundary value expansions”. 
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Although the classical Lie point symmetry classification of Q, with D[6) 
and R{0) arbitrary, selects only constant or unbounded reaction functions, the 
complete nonclassical symmetry classification of reaction-diffusion equations, 
admits a broader range of possibilities [iiiniiiiiiii], even admitting the Arrhe¬ 
nius reaction term |12j . Nonclassical symmetries, in the sense of Bluman and 
Cole |13j . leave invariant a system consisting of the original governing equation 
Q plus the invariant surface condition that restricts the solution set to only in¬ 
variant solutions. The concept of nonclassical symmetry extends more generally 
to that of compatibility with an invariance condition, that may be of higher or¬ 
der than a point transformation or a contact transformation El- Nonclassical 
symmetry analysis reveals that with a suitable nonlinear diffusivity function, 
after a change of variables, the Arrhenius reaction-diffusion equation admits 
separation of variables, resulting in a linear system. In the following Section 2, 
time-dependent radial solutions of this system are constructed in two and three 
dimensions. In the first example, the complicated nonlinear diffusivity function, 
as well as the temperature, are given exactly. With the same strongly increas¬ 
ing diffusivity, the exponentially heating solution is mathematically equivalent 
to that involving conduction through a finite domain with an endothermic re¬ 
action and exponential cooling. The strongly increasing diffusivity resembles 
soil-water diffusivity, with the sink term representing plant root absorption. In 
the case of ideal maximal cooling at the boundary, the nonlinear diffusivity is 
bounded and it is the fixed point of a rapidly converging contraction map. The 
similarity form of the Kirchhoff variable is given exactly and it is asymptotically 
of the same form as the temperature. The diffusivity varies so little that the 
exact solution for the Kirchhoff variable closely approximates the temperature 
at all times. Particular attention is paid to the case of Newton cooling at the 
boundary, with non-zero Biot number. 

In Section 3, it is shown that the same type of solution construction applies to 
a class of nonlinear reaction-diffusion equations with spatially varying diffusivity. 
The results of Sections 2 are extended to allow for spatially variable diffusivity 
and a nonlinear sink term that applies to water transport in unsaturated soil 
with extraction by plant roots. 


2 Nonclassical reduction of nonlinear reaction- 
diffusion to the Helmholtz equation 

A full nonclassical symmetry classification of nonlinear diffusion-reaction equa¬ 
tions in two spatial dimensions, was given in m- Firstly, Q is expressed in 
terms of the Kirchhoff variable (e.g. |I5)L 

u = uq + f D{9)d9 . (3) 

Jo 
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If Mo = — Jq° D{9)d9 for some 9o > 0, then 


M = / D{9)d9 . 
ho 

In that case, a boundary condition u = 0 corresponds to 9 — 9o. 

In terms of the Kirchhoff variable, the reaction-diffusion equation is 

F)i I 

F(m)^=V2m + Q(m), 

where Q{u) = R{9) and F{u) = 1/D{9). The starting point of the reduction 
to the Helmholtz equation is the observation that Equation Q has a simple 
nonclassical symmetry 


( 4 ) 


M = ; t = t + e ; 

whenever F and Q are related by 

Q(u) = AuF{u) -\- KM, 


with e € 


(5) 

( 6 ) 


for some H, k G K. This is not a classical symmetry because it does not in 
general leave the equation Q invariant, except when one also assumes the 
invariant surface condition ut = Au. Then the symmetry reduction leads to 

= 0 with M = e'^‘$(x). (7) 


In terms of the original temperature variable 9, the relation (§ is 


R{9) 


r A 1 

r /•'--! 

ml 

uq-\- dd 

Jo 


( 8 ) 


This gives an explicit construction of R{9) from D{9). Some basic combinations 
{D{9)^R{9)), with i?(0) = 0, are given in Table 1. Note that the R{9) function 
in case (d) of Table 1 agrees asymptotically with the Arrhenius reaction term 
as inverse temperature approaches oo. 

However, in order to construct D{9) from R{9), it is necessary to solve a 
first-order nonlinear differential equation. From §, 


D{9) = u’{9) 


Au 

R{9) — Ku 


(9) 


By making u{9) the subject, then differentiating, this implies 


AD'(9) = 


m 


R 


DA 


A-R'{9) 

R ■ 


( 10 ) 


The equation 0 for $ is simply the linear Helmholtz equation if k = > 0, 
and the Laplace equation if k = 0. If k = —K^ < 0 with R > 0, from ([^ A 
must be positive. Then the equation for $ is the modihed Helmholtz equation 
for which the radial solutions are unbounded either at the origin r = 0 or at 
infinity. 
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D{ 0 ) 

m 

(a) 

0 "^ (to > -I) 

« Qm+l ^ 0 

TO -f 1 TO -f 1 

(b) 

e® 

K[e« - 1] - A[e-^ - I] 

(c) 

cosh(0) 

Ksinh(0) — Atanh(0) 

(d) 


Bo^(B + 6»)e-'^/®(Roe"'®/® - AB) 

kB \ 0 ) K 

RnB(B + 0)e-^/^ -AB20 


Table 1: Reaction-diffusion combinations that admit nonclassical scaling sym¬ 
metry 

2.1 An explicit solution with Arrhenius reaction: k = 0 

In the case k = 0, equation is linear homogeneous, allowing direct integration 
to obtain 

-1 

with 

C(«) = j|^exp(/ (12) 

where ci is an arbitrary constant that is positive (negative) if i? is a strictly 
positive (negative) source (sink) term. However the spatial dependence of it, 
given in 01 is then governed by Laplace’s equation. For example, the isotropic 
radial solutions it(x, t) = u{r, t) in two or three dimensions can only be a con¬ 
stant added to the singular point-source(sink) solutions with the source(sink) 
strength varying exponentially in time. The isotherms are specified exactly by 
the mapping that follows from Equation ( |11[ ): 

(t, 6*) u !->■ $ = e~^*u !->■ r. 

In two and three dimensions, the radial solutions take the form r = exp( ) 
and r = respectively. 

For the case of the Arrhenius reaction Q, the solution 0 is valid when the 
diffusivity is exactly 

D = exp(R/6») exp ("^6»exp(R/6») - ^E,{B/0) \ . 

where Ei is the exponential integral. When 6 is large, we can use the power 
series for the exponential integral 

oo ^ 

E^{x) = 7 -f ln(a;) + X! ^ ’ 
nn\ 
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where 7 = Euler’s constant, to find 

JJ ^ /Rq j^-AB/R q qAB/Rq ^Ae/R q 

Rq 

When 0 is small, we can use the asymptotic expression 


E,{x) ~ - 

X 


n—1 


1 +V — 

/ ^ j,n 


to find that Z? —>■ 0 as 0 —>■ 0 . 

This nonlinear diffusivity is plotted in Figure lb alongside the Arrhenius 
reaction function, plotted in Figure la. 




Figure 1: (a) Exact temperature-dependent diffusivity allowing separation of 
variables of Arrhenius reaction-diffusion, (b) Arrhenius reaction rate, as a func¬ 
tion of temperature. 

Example 1. A heat conduction problem with Arrhenius endothermic 
reaction term 

Even though the enthalpy of reaction is negative, an endothermic reaction 
will proceed naturally if the reaction products result in a lowering of the Gibbs 
free energy. An endothermic reaction may have a single activation energy, so 
that the rate of heat absorption is described by an Arrhenius function. Consider 
a spherical vessel containing the reagents of an endothermic reaction. Rq is 
negative. A small spherical decaying radioactive heat source at the centre, 
supplies a quantity of heat Q. Then, since u > 0 and D > 0 when 0 > 0, 
now imply ci < 0 and A < 0. A solution 

u = e-l^l‘[c2-^] 

r 

can be made to satisfy the boundary conditions 

u = 0 at r = ri, 

—inr^Ur = at r = rg. 

Note that the latter is simply a heat flux condition that can always be expressed 
as a linear condition in terms of the Kirchhoff variable u. 
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Example 2. Porous media flow with plant-root sink term 

The combination of exponentially growing nonlinear diffusivity, with bounded 
sink term, applies to water transport through a web of plant roots [16]. The 
domain of the problem is considered to be a cylindrical pile of soil of radius 
r = ri with a vertical cylindrical injection well in the centre (r = rp < ri). The 
dependent variable 9 now designates the water content above the plants’ wilting 
point, where the sink term approaches zero as roots fail to draw water. In this 
application, Rq is negative. Then, since u > 0 and D > 0 when 6 > 0 , 
now imply ci < 0 and A < 0. A solution 

u = — C 3 Inr] 

can be made to satisfy the boundary conditions 

u = 0 at r = ri, 

—2TTrUr = at r = rp. 

These boundary conditions allow for injection of total water volume Q per unit 
length of injection well into a large cylindrical soil mound that is exposed to the 
soil-controlled second stage of evaporation (e.g. m) at its outer boundary. 

This solution applies analogously to a vertical current-carrying wire along 
the axis of a cylindrical region, supplying energy for an endothermic reaction. 

2.2 Construction of bounded temperature solutions: k > 0 

Example 3. Heat conduction with exothermic reaction in a compact 
region 

With K = > 0, u = d)(r) exp(At) satisfies the linear Helmholtz equation. 

A non-negative isotropic solution can satisfy boundary conditions 

Ur = 0 , r = 0 , 
u = 0, r = ri, 

by choosing $ = jo(Arr) in three dimensions, $ = Jo{Kr) in two dimensions 
and $ = cos{Kr) in one spatial dimension, with K = Ai/ri, Ai being the first 
zero of the spherical Bessel function jg (-^i = '^)i the standard Bessel function 
Jg (Ai « 2.4048) or the cosine function (Ai = 7 r/ 2 ) respectively. Further, 
from the separation of variables in Q , that same solution satisfies other linear 
homogeneous boundary conditions such as 

—Ur = Bi u, r = r 2 < ri, 

with Bi constant. This resembles the linear condition for Newton cooling, in 
which Bi is the Biot number, after rescaling and non-dimensionalising vari¬ 
ables. The left hand side is heat flux but the right hand side is a multiple 
of the Kirchhoff variable, rather than temperature. Below, it is demonstrated 
that the nonlinear diffusivity has bounded variation and is close to being con¬ 
stant. This means that the above boundary condition is very close to that of 
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Newtonian cooling to a very low-temperature environment. If the physical pa¬ 
rameters r 2 , D{0) and Bi are prescribed, then the solution parameters are given 
\yy A = —K^D{0), where K is the unique solution of —^'{Kr 2 )/^ — Bi/K. K 
must lie between 0, where —^'{Kr 2 )/^ = 0 and Ai/r 2 , where —^'{Kr 2 )/^ ap¬ 
proaches infinity. With these homogeneous boundary conditions, there always 
exists an extinguishing solution in which the temperature approaches zero uni¬ 
formly everywhere. A demonstration of the stability of this similarity solution 
is given in the Appendix A. 

In the case k > 0, is equivalent to the canonical form for an Abel equation 
of the second kind, via 


w = Ku — R{9), z = —A6 — R{9), 
dw AR{9) 


(14) 


The standard list of known integrable forms of the Abel equation is given in 
[TS] . In principle, this list may be used to furnish R{9) functions that produce 
solvable forms of equation ( |l4| ). If g{z) a particular solution for an integrable 
form of (14), and g~^{w) is the corresponding inverse function, an R(9) function 


that produces (14) for w{z) is given implicitly according to 


9 = 


R + g-\-R) 
-A 


(15) 


It is not clear if the segmented solution method of m is practicable for the 
general case. 


A good approximate analytic reconstruction of D{9) may be obtained by 
applying a contraction map towards a fixed point. With uq = 0, the solution to 
satisfying initial value D{9) —> 0, 0 —>■ 0, must be a fixed point of the map 


Dr^+l{9) = MDr,{9) 


-A Dn{s)ds 
D^{s)ds - R{9)' 


When considering the separated solution Q, |A| and K may be set to 1 by 
choosing dimensionless length and time coordinates Kr and \A\t. 


Dn+l{9) = 


Dn{9) 


D^i9)-Ri9)/9^ 
where Dn{9) is the running mean value. 


(16) 


1 r 

Dn{9) = -J D^{s)ds. (17) 

In the case of the Arrhenius reaction term, we may set R{9) to be Roe~^^^ 
by using a dimensionless temperature variable 9/B. Then the maximum value 










of R{0)/9 is Ro/e. 


From (16), using D{0) = 11(0), it follows that D{0) = 1. From (16) it also 
follows that D{9) > 1 for 0 > 0. If D{9) is bounded for 9 € (0, oo) and mono¬ 
tonic for 9 sufficiently large, D(9) must have a limit D{oo) = D{oo). In that 
case it follows from (16) that D{oo) = 1. 


Consider AI as a mapping on the set S of bounded continuous functions / 
on (0, oo) such that / — 1 is non-negative, 

S = {f G (71(0, oo) : 36 Vcc e (0, oo) 0 < /(x) - 1 < 6}. 


Suppose that Dq,Eq G {5}, = M'^Dq, E„ = M'^Eq, where if„+i(6) and 

En{9) are defined in a analogous way to (16) and 0 above. Now 


|I3„+i(6)-F;„+i(6)| 


Dn{x)dx /g En{x)dx 

/q Dn{x)dx - R{9) En{x)dx - R(9) 

9Dn _ 9En 

9Dn - R{9) 9En - R{9) 

(Dn - Er,)R{9)/9 
[D^- R{9)/9][E^- R{9)/9] ' 


If i?o < e, then inf{Zl(6)} = 1 > sup{i?(6)/6} = i?o/e, consequently 

\\En+i - Du+i\\oo ^ _ sup{i?(6)/0} _ 

\\En-Dn\\ao ~ [inf{Il„(6»)} - sup{i?(6»)/6»}][inf{F;„(6»)} - sup{i?(6»)/6»}] ■ 

||.^n-t-l IIqq ^ 1 

\\En — Dn\\ao tjR q — 1 Rf^ je 

This is a contraction map, provided 


i?o < 


(3 - y5)e 
2 


1.0383 . 


For example, consider the case Rq = 1 for which the above estimate of the 
contraction factor is 1/1.086. In practice, the convergence rate is better than this 
modest value suggests. In the case R{9) = exp(—1/6), since D{0) = D{oo) = 1, 
it is natural to choose Dq{9) = 1. Using the standard notation /? = 1/6, the 
iterative estimates are 


Do 

Di 

D2 

I 


1 , 

- OO 

(noting I3e-^ < e"! < 1), 

^ n=0 

I 


/-e-t 

“ Hi 


L 


dl3 = p- 


in-2)\ n^/3 




/c =0 


kpk 

k\ 
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Figure 2: D{9) constructed by series about 9 = 0,0.5,10, as well as approxima¬ 
tions Hi ( 6 »), H 2 ( 6 I) 


Each term of the I summation is bounded above by (n — 2)!/n"“^, which is of 
order e~'^ jas n ^ oo. The n = 5, terms have a combined value of less 
than 0 . 01 . 

Note that the sum of the terms has a maximum value of less than 0.01. 
In Appendix B, D{9) is calculated accurately after obtaining exact series forms 
for u{9). The simpler approximate diffusivity functions Dj{9) are shown in 
Figure to agree well with the exact solution. 

The diffusivity function has a single local maximum where its value is no 
higher than 45% above its mean value of 1. From (10), a stationary point of 
D{9) may occur only where D{9) = 1 + R'{9) = l + RQ9~‘^e~^^^ . This expression 
has a maximum value, implying 


Dm < 1 -f dAo/e^ « 1.5413 


(18) 


In effect, since at small-t, local disturbances extend by diffusion to a depth 
proportional to and varies by only 20 %, the diffusivity is effectively 
constant for some practical purposes. This is seen in the solution, plotted in 
Figure 3, wherein the temperature 9{r, t) asymptotically approaches the Kirch- 
hoff variable u{r,t) which is identical to the temperature when D = 1. 

The construction of the nonlinear diffusivity function by a contraction map 
relied on the fact that with Arrhenius reaction, /3Q(/1) (= R{9)/9) has an upper 
bound less than 1. The same construction will apply to other reaction laws with 
this property, for example Q{/3) = /l'"exp(—/I) with m > —1. This includes the 
case m = — 1/2 that follows from the kinetic theory of gases. 
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Figure 3: Temperature profiles at times |^|t=-1.5,0,l-5 and 2.5 (solid) approach¬ 
ing the Bessel function profile of u(r,t) (x crosses) from below, as time increases. 

2.3 The case k < 0. 

If K = —K^ < 0, from § with R{9) > 0, A must be positive. The diffusivity, 
D{6), has a singularity where at leading order D'{9) ^ {k^/A)D^ and D ~ 
{&a — 9)~^^‘^ at some freely chosen positive value 9 = 9a. In the modelling 
of unsaturated soil-water flow, the location of the singularity is chosen to be 
slightly greater than the volumetric water content at saturation m- 

Also the equation for $ is the modified Helmholtz equation for which the 
radial solutions are unbounded either at the origin r = 0 or at infinity. There 
is an exact solution for exothermic reaction-diffusion on the domain r > tq, 

u = e^^Ko{Kr), (19) 

where Kq is the modified Bessel function of the second kind. Although u is 
unbounded, the corresponding value of 9 remains less than the singular value 
0a, due to the singularity in D{9) within the integrand in the Kirchhoff trans¬ 
formation 9 ^ u. This singular nonlinear diffusivity does not occur in heat 
conduction but it may occur in models of unsaturated soil-water flow. However, 
in that application, positive distributed water sources are not common. It is 
more common to consider distributed sinks, due to plant roots. 

Example 4. Porous media flow with plant-root sink term: k < 0 

As in Example 2, we consider a distributed sink with Rq < 0, due to plant 
roots. Because of plant root extraction, A < 0 so that a finite amount of wa¬ 
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ter is supplied over all time through the inner surface at r = tq. From Q, 
D > 0, A < 0 and k < 0 implies R < 0. However in this case, from the solution 
(19), 0 does not reach zero at any finite value of r but instead approaches zero 
at infinite distance. It is to be noted that the modified Helmholtz equation has 
previously been applied to steady state solutions of unsaturated soil-water flow, 
by analogy to external problems of Helmholtz acoustic scattering Just as 
in the latter, we may construct analogous exterior solutions for boundaries of 
various shapes. However, at large values of r it is well known that the solu¬ 
tions agree asymptotically with the isotropic ones, so the spherical or circular 
boundaries are canonical. 


3 Incorporating heterogeneity 

We now consider the case in which the medium (or substrate) is spatially hetero¬ 
geneous. The appropriate governing reaction-diffusion equation can be written 


0t = V.[/(x)Z?(0)V0]+i?(0), (20) 


where x S fl C K". Here, the heterogeneity is represented by a positive differen¬ 
tiable function / that is the amplification factor of the diffusivity. By rewriting 
this equation in terms of the Kirchhoff variable (|^, we obtain 

F(u)ut = V.[/(x)Vu] -I- Q{u), 


(the heterogeneous equation corresponding to (|^) where Q{v) = R{9) and 
F{u) = 1/D{6), as before. Equation (20) admits the same simple nonclassical 
symmetry ([^, whenever F(u) and Q{u) are related by Equation ®> .This means 
that R{9) can be constructed from D{9) using ([^, or D{9) could be constructed 


from R{9) by solving 0- By again setting m(x, t) = 
second order linear equation for $(x). 




$(x) we derive the 


/(x)V2$(x) + V/(x).V$(x) -b K$(x) = 0, (21) 


which is valid in n dimensions. In the case where A < 0, the solution for 9 
will be asymptotically similar to u. The relationship between D{9) and R{9) is 
the same as that described in the preceeding section, and the solutions for D{9) 
found about also hold in the hetergeneous examples described below. 


3.1 The case /t = 0 

Example 5. Porous media flow with a heterogeneous substrate: k = 0 

As an example, we now reconsider the porous medium plant-root extraction 
problem described in Section 2.1 above, with the inclusion of spatial heterogene¬ 
ity representing a changing scale of soil pores. In Miller self-similar soils, / is 
the geometric scale factor of the soil pores ([23l|24]), with /(ro) = 1 at some 
chosen reference point. 


12 




In this example, the problem is written in (one-dimensional) radially sym¬ 
metric coordinates and the heterogeneity is described by /(x) = f{r). The 
differential equation for <i)(r) replacing $(x) in (21), is 


r dr 


dr 


dr dr ’ 


which can be directly integrated for arbitrary /(r) to give 

1 


u{r,t) = e^*^{r) = ce 




rf{r\ 


-dr 


with c constant. From ([T^), we deduce that A < 0. Having freely chosen /(ro) 
to be 1, the solution, in terms of the Kirchhoff variable, satisfying boundary 
conditions (13), is 


u{r t) = —^ _ — 

’ 27r X sfis) 


ds. 


3.2 The case /t > 0 


Example 6. Heat conduction in a heterogeneous medium 

Let K = K^. The differential equation for <I>(r) replacing $(x) in (21) is 


f d f d^\ dfd<^> ^ ^ 

rTr 


( 22 ) 


With /(r) any power-law of the cylindrical radius, exact solutions are readily 
available (e.g. m)- 

With heterogeneity described by the factor /(r) = ro/r, (22) reduces to 
Airy’s equation, with general solution 


When /(r) is of the form 

( 23 ) 

the resulting ODE is merely a homogeneous Euler equation, with general solu¬ 
tion 

$(r) = -I- . 

This solution can be made to satisfy boundary condition <i)(ri) = 0 by setting 

2y'l-{Kroy 
C2 = -Ciri'' 

For the case with Kr^ > 1, the general real solution is 

$(r) = - [ci cos(wlogr) -|- C 2 sin(a;logr)]; uj = \/ {Kro)'^ — 1. 

For the above cases, solution parameters may be chosen so that u(ri) = 0. 
However, r = 0 is a singular point where thermal conductivity is either zero or 
infinite. The solution may be regarded as exterior to the surface of a hot wire 
at r = rg from where a finite quantity of heat is supplied. 
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3.3 The case k < 0 


Example 7. Porous media flow with a heterogeneous substrate: k < 0 

With the same heterogeneity as in the previous model (23) but with k = 
—K^ < 0, the solution with /(r) = (r/ro)^ is simply 


$(r) = cir“^+v'i+(^’'o)= c2r“^“Vi+(^’'o)2 


and with /(r) = ro/r, 

The parameters q and K may be chosen so that the exterior solution satisfies the 
same boundary conditions on [ro,oo] as in the case of a homogeneous medium. 


4 Conclusion 

Provided the nonlinear diffusivity and the nonlinear reaction term satisfy a 
single relationship, the Kirchhoff variable u, which is the integral of nonlinear 
diffusivity, admits solutions that are obtainable by separation of variables to a 
linear system, whose solution is an exponential in time, multiplying an arbitrary 
solution of the Helmholtz, modified Helmholtz or Laplace equation in space. The 
heat flux is merely — Vu, which is given explicitly everywhere in these solutions. 

If the nonlinear diffusivity function is specified, then the compatible reaction 
function can be constructed directly by integration. If the nonlinear reaction 
function is specified, then the diffusivity function is the solution of a differen¬ 
tial equation that is equivalent to an Abel equation if u satisfies a Helmholtz 
equation, or to a separable equation if u satisfies Laplace’s equation. To the 
best of our knowledge, this provides the only known exact closed-form solutions 
in two and three dimensions, of nonlinear reaction-diffusion equations with the 
classical Arrhenius reaction term. Even when D(0) satisfies an Abel equation 
that is not of known solvable type, in some circumstances it may be specified 
to arbitrary precision using exact series expansions. 

This construction is valid in any natural number n of dimensions, and it 
generalises also to heterogeneous extensions of the Helmholtz factor equation. 
Applications are given for radial solutions of exothermic reactions with nonlinear 
heat conduction, endothermic reactions with nonlinear heat conduction, and 
water flow from a supply well into a cylindrical soil mound with soil-limited 
evaporation at the outer boundary. The logistic shape of the negative Arrhenius 
reaction term resembles the behaviour of distributed plant roots that have a 
maximum and minimum value of water uptake rate near saturation and wilting 
point respectively. 

As is well known from acoustic scattering theory, exterior solutions of the 
Helmholtz equation typically asymptotically approach radial symmetric solu¬ 
tions at large distances from the scattering surface. Hence, the radial solutions 
illustrated here are in a sense, canonical. The solution method used here, in¬ 
volves a free function of the Helmholtz equation so it would not be difficult to 
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use known non-radial solutions. The special solutions presented here could at 
least be used as bench tests for more broadly applicable approximate solution 
methods. 

This approach will lead to ongoing investigations of other nonlinear partial 
differential equations with more than one free function, that may admit special 
nonclassical symmetry reductions. 
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Appendix A. Stability of similarity solution. 

It is shown here that the extinguishing radial similarity solution of the Arrhe¬ 
nius reaction-diffusion equation given in Section 2.2, is exponentially stable to 
small perturbations. Change variables to a set of canonical variables of the 
nonclassical symmetry, namely r,t and v = ue~^*. In terms of these variables, 
the reaction-diffusion equation may be expressed 

F{u)vt = -I- K^v, 

with boundary condition v{ri) = 0. The similarity solution is a pseudo-steady 
state. Vs = 4’(r) = Aq Jo(Air/ri), corresponding to exponential decrease of the 
Kirchhoff variable Us = and satisfying $'(0) = 0 and $(ri) = 0. Now 

consider a perturbed solution, in plane polar coordinates 

V = <i)(r) -I- w{r, 4>, t), 

with twice-differentiable initial condition 

V = $(r) -I- wo(r, (j))-, [wol < e << 1- 

Then order-e perturbation w satisfies the same homogeneous boundary condi¬ 
tions, plus 27r-periodicity with respect to (j), plus 

F (us(r, t)) Wt = K^w. 

Unlike the original equation for temperature 9, this is a linear equation with no 
squared derivative terms, allowing recourse to comparison theory of reaction- 
diffusion equations (e.g. [IS]). Note that 1/F{u) = D{0) and that from Section 
2.2, Dm > D > D{0) = —A/K"^. Hence by comparison, \w{r,t)\ must decay at 
least as fast as the solution to the linear initial-boundary problem with smaller 
diffusivity and larger reaction term, 

Qt = DiO)V^Q + DmK^Q ; Qir,(f,0)=wo{r,(fy, Q(ri, </., t) = 0. 

Let 

T = \A\t/K‘^; and q = 

Then we have the standard linear heat equation with q satisfying 

the same initial and boundary conditions as Q. The solution g(r, (j), t) may 
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be expanded as a standard Fourier-Bessel series (e.g. [IS]). Without loss of 
generality, we neglect the first component in the series expansion of wq , which is 
a multiple of and which takes the form of the assumed unperturbed 

similarity solution for v. Each of the other terms in the series for Q is of the 
form 

Jn{^n,-mr/ri) cos(n^) + Bn,m sin(n^)] ^ ( 24 ) 


where n S N, Xn,m is the m’th zero of Bessel function Jn, while An^m and Bn^m 
are arbitrary real coefficients. Note that beyond Ao,i, the next smallest root is 
Ai.i. From the estimate (18), each of the terms (24) is decreasing exponentially 
in time, provided 


Ro < e2[(Ai.i/Ao.i)2 - l]/4 « 2.8425. 


(25) 


In terms of the original physical parameters of the unsealed boundary value 
problem, this criterion is 


-^0^1 

BXI,D{0) 


< 2.8425, 


(26) 


which is satisfied in practical cases. 


Appendix B. Evaluating D{9) using asymptotic series. 


The reaction term R{0) = Rq exp(—5/0) is of special interest. With k = K^, 
the parameters K, |A| and B in equation may be set to 1 by adopting 
appropriate dimensionless variables. By successive isolation of leading-order 
terms, the asymptotic behaviour of u{0) near 9 — 0 can be shown to be 


m( 0) ~ —sign(A)0-I-i?o0exp(—1/0) ^ 0 exp(—r/0) ^ 0'"gr,m- (27) 

r—0 m—0 


This series structure assumes mq = 0 but extension to mq > 0 should be straight¬ 
forward. Adopting Qr^-i = 0, substitution into © shows that 


<?o.m = (-l)™m!, for m > 0 ; 
(-sign(A)i?o)’ 


Qr,0 — 

0 r+l.m-|-l = I 


r -I- 1 

sign(A)i?o f sign(A)(r — m 


for r > 0 ; 

0r-t-l,m T (r j l)0r,r; 


5 , 


0 

r m 


- {r + l)qr,m+l + X! X! 9r-s,m-i [(s -f l)qs,l “ (s “ l)qs,i-i] 
for r, TO > 0 . 


s =0 1=0 


(28) 

(29) 


(30) 


With the set of coefficients {qr,m} determined iteratively as above, the series 
(27) is divergent, but can still be used to give accurate specifications of u(0) for 
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6 sufficiently small. Partial sums of (271 may be produced by evaluating both 


r and m sums up to some maximum integer. Performance of the partial sums 
is improved by substituting the known asymptotic form 




(31) 


771=0 


A Taylor series expansion for u{9) centered on a non-singular point 0 < 0o < oo 
is comparatively easy to derive. We adopt 


i(y) = ^"2/” 


y = 


n—0 


00 - 0 , 
00 ’ 


(32) 


and by substitution, can show that 

1 


Ai+i — 


{i + l)(Ao - i?oexp(-l/uo)) ( 


i sign(A) 6 »oA, 


— — n + l)Ai_„+i A„ -f -^ 1 ^ 1(71 + 1, 2; —I/Oq) 

for z > 0 . 


(33) 


Here Aq = u{9o) is assumed known, and the above sum evaluates to zero for 
z = 0. The asymptotic series form of u{9) as 0 —>■ 00 can also be ascertained 
and appears to have a non-zero radius of convergence. 

We can use the above series to evaluate zz for 0 < 0 < 0max to a prescribed 
accuracy by matching expansions about different points. To demonstrate, take 
i?o = 1 and A < 0. Using (27) with rmax = nT-max = 10, we can ascertain 
u(l/10) = 0.1000041579094 to 12 significant figures. A Taylor series expansion 
about 9q = 1/2, with u{l/2) =0.55582409195937 then covers the domain 0.1 < 
9 < 0.9, and matches the accuracy of the known value at 6* = 1/10. Introducing 
a second expansion about 9o = 10, with zz(lO) = 11.79313028084656, extends 
the domain of the solution to 9 = 19.5, and is compatible with the known value 
u(9/10) = 1.1135172087801. Figure 2 shows the resulting diffusivity D{9) = 
du/d9. To obtain a solution of greater accuracy, we can begin with evaluation 


at a point 9 < 1/10 according to series (27). A greater number of linked Taylor 


series may then be needed to cover 0 < 0* < 19.5 as done above. 
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